clear; clc;
import Algorithm.BlochEquation.ParametricBlochEquation

%%
T1=inf; T2=10e-6;
model = ParametricBlochEquation(T1, T2, ...
                           'ratio_x', 3.0, ... % ratio of `control amplitude' to `control frequency'
                           'ratio_kc', 1.2, ...% ratio of `control frequency' to `Lamor frquency'
                           'pumping_rate', 10e3);

model.set_Omega('omegaX', 2*pi* 0.0,  ... 
                'omegaY', 2*pi* 10.0, ...
                'omegaZ', 2*pi* 34e3);
        
sol_ana = model.solveAnalytically;
sol_num = model.solveNumerically;

plot(sol_ana.mx, 'rd', sol_ana.my, 'bo', ...
     sol_num.mx, 'r-', sol_num.my, 'b-');
legend({'mx', 'my'});